Wstęp do uczenia maszynowego

Projekt nr 2

Jakub Piwko, Malwina Wojewoda

Pakiety

#!pip install mca
#!pip install prince
import pandas as pd
import numpy as np
#from sklearn.cluster import KMeans, AgglomerativeClustering
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import MinMaxScaler
from sklearn.decomposition import PCA
import mca
import prince
import itertools
import scipy.stats as ss

Wczytanie danych oraz ich opis

Jak dowiadujemy się z opisu danych są one obrobioną nieco wersją zbioru USCensus1990raw. Zawierał on 1 % próbki z pełnego spisu powszechnego z 1990r. W naszym pliku rekordy zostały losowo przemieszane. Aby uzyskać zbiór, którym się zajmujemy zostały wykonane następujące operacje:

• randomizacja: rekordy z USCensus1990raw zostały losowo permutowane
• zmiany w kolumnach: kilka zmiennych ciągłych zostało zdyskretyzowane, a zmienne dyskretne, które miały dużo wartości zostało zwinięte do mniejszej liczby wartości; kolumny, które nie uległy zmianie oznaczone są prefiksem "i", a te które były zmienione - prefiksem "d"; dokładne zmieny zostały opisane w pliku USCensus1990.mapping.sql
• usunięto kolumny, które były mało przydatne

Można więc powiedzieć, że na danych został już wykonany preprocessing.

Pierwszym atrybutem jest caseid i powinnien on być ignorowany podczas analizy.
Wszystkie dane są numeryczne.

X = pd.read_csv('USCensus1990.data.txt')
X.head()

	caseid	dAge	dAncstry1	dAncstry2	iAvail	iCitizen	iClass	dDepart	iDisabl1	iDisabl2	...	iTmpabsnt	dTravtime	iVietnam	dWeek89	iWork89	iWorklwk	iWWII	iYearsch	iYearwrk	dYrsserv
0	10000	5	0	1	0	0	5	3	2	2	...	0	5	0	2	1	1	0	11	1	0
1	10001	6	1	1	0	0	7	5	2	2	...	0	1	0	2	1	1	0	5	1	0
2	10002	3	1	2	0	0	7	4	2	2	...	0	2	0	2	1	1	0	10	1	0
3	10003	4	1	2	0	0	1	3	2	2	...	0	1	0	1	1	1	0	10	1	0
4	10004	7	1	1	0	0	0	0	2	2	...	3	0	0	0	2	2	0	5	6	0

5 rows × 69 columns

X.shape

(2458285, 69)

Sprawdźmy czy w zbiorze są jakieś braki danych:

X.isna().sum().sum()

0

X.isnull().sum().sum()

0

Wygląda na to, że w naszym zbiorze nie ma braków danych.

X.duplicated().sum()

0

Nie ma też zduplikowanych wierszy.

Sprawdzimy jeszcze jakie są unikalne wartości w każdej z kolumn:

pd.set_option('display.max_rows', None) 
print(X.apply(lambda col: col.unique()))

caseid       [10000, 10001, 10002, 10003, 10004, 10005, 100...
dAge                                  [5, 6, 3, 4, 7, 1, 2, 0]
dAncstry1               [0, 1, 11, 3, 2, 9, 8, 10, 4, 6, 7, 5]
dAncstry2              [1, 2, 12, 3, 4, 10, 7, 5, 11, 8, 6, 9]
iAvail                                         [0, 3, 4, 2, 1]
iCitizen                                       [0, 3, 4, 2, 1]
iClass                          [5, 7, 1, 0, 6, 3, 4, 2, 8, 9]
dDepart                                     [3, 5, 4, 0, 2, 1]
iDisabl1                                             [2, 0, 1]
iDisabl2                                             [2, 0, 1]
iEnglish                                       [1, 0, 2, 3, 4]
iFeb55                                                  [0, 1]
iFertil         [1, 3, 0, 4, 7, 2, 5, 6, 12, 8, 13, 10, 9, 11]
dHispanic                       [0, 1, 9, 2, 6, 3, 4, 8, 7, 5]
dHour89                                     [4, 1, 3, 0, 5, 2]
dHours                                      [3, 1, 4, 0, 5, 2]
iImmigr                     [0, 5, 1, 8, 3, 10, 6, 7, 4, 2, 9]
dIncome1                                       [2, 1, 0, 3, 4]
dIncome2                                                [0, 1]
dIncome3                                                [0, 1]
dIncome4                                                [1, 0]
dIncome5                                                [0, 1]
dIncome6                                                [0, 1]
dIncome7                                                [0, 1]
dIncome8                                                [0, 1]
dIndustry           [10, 4, 1, 0, 9, 7, 6, 3, 5, 8, 11, 12, 2]
iKorean                                                 [0, 1]
iLang1                                               [1, 2, 0]
iLooking                                             [0, 2, 1]
iMarital                                       [1, 0, 4, 2, 3]
iMay75880                                               [0, 1]
iMeans              [1, 10, 0, 11, 2, 4, 8, 5, 12, 3, 7, 9, 6]
iMilitary                                      [4, 0, 2, 3, 1]
iMobility                                            [2, 1, 0]
iMobillim                                            [2, 0, 1]
dOccup                             [3, 2, 4, 0, 6, 1, 5, 8, 7]
iOthrserv                                               [0, 1]
iPerscare                                            [2, 0, 1]
dPOB                                     [0, 3, 4, 6, 2, 1, 5]
dPoverty                                             [2, 1, 0]
dPwgt1                                            [1, 2, 0, 3]
iRagechld                                      [4, 2, 0, 3, 1]
dRearning                                   [3, 2, 0, 4, 5, 1]
iRelat1         [0, 1, 2, 10, 12, 9, 7, 4, 3, 13, 11, 6, 8, 5]
iRelat2                         [0, 6, 9, 1, 8, 2, 4, 3, 7, 5]
iRemplpar    [0, 121, 112, 211, 111, 122, 222, 221, 223, 11...
iRiders                            [3, 1, 0, 2, 4, 8, 6, 5, 7]
iRlabor                                  [1, 6, 0, 3, 4, 2, 5]
iRownchld                                               [0, 1]
dRpincome                                   [3, 2, 0, 4, 5, 1]
iRPOB        [22, 10, 21, 23, 51, 24, 52, 40, 31, 35, 36, 3...
iRrelchld                                               [0, 1]
iRspouse                                 [3, 1, 6, 4, 0, 5, 2]
iRvetserv               [0, 1, 3, 8, 10, 4, 7, 2, 5, 9, 6, 11]
iSchool                                           [1, 2, 3, 0]
iSept80                                                 [0, 1]
iSex                                                    [1, 0]
iSubfam1                                          [0, 2, 3, 1]
iSubfam2                                          [0, 1, 2, 3]
iTmpabsnt                                         [0, 3, 1, 2]
dTravtime                                [5, 1, 2, 0, 4, 3, 6]
iVietnam                                                [0, 1]
dWeek89                                              [2, 1, 0]
iWork89                                              [1, 2, 0]
iWorklwk                                             [1, 2, 0]
iWWII                                                   [0, 1]
iYearsch     [11, 5, 10, 4, 8, 1, 7, 6, 12, 14, 15, 2, 0, 9...
iYearwrk                              [1, 6, 0, 7, 4, 2, 3, 5]
dYrsserv                                             [0, 1, 2]
dtype: object

pd.set_option('display.max_rows', 10) #back to default

Oryginalny zbiór jest duży, ograniczymy go do losowej jednoprocentowej próbki, dodatkowo usuniemy zmienną caseid.

df = X.sample(frac=0.01, random_state=17)
df.drop('caseid', inplace = True, axis = 1)

df.shape

(24583, 68)

Analiza zmiennych

Przejdziemy teraz do głębszej analizy zmiennych zawartych w naszej ramce. Przyjrzymy sie mapie korelacji a także zależnościom między zmiennymi.

plt.figure(figsize=(40, 40))
sns.heatmap(df.corr(method = 'pearson'), annot=True)
plt.show()

Z powyższej mapy możemy odczytać, że w ramce danych występują wysoko skorelowane zmienne. Podzbędziemy się tych o najwyższym module wartości korelacji (czyli większym od 0.7).

correlated_features = set()
corr = df.corr(method = 'spearman')
for i in range(len(corr .columns)):
    for j in range(i):
        if abs(corr.iloc[i, j]) >= 0.70:
            colname = corr.columns[i]
            correlated_features.add(colname)

print(len(correlated_features))
print(correlated_features)

29
{'iMeans', 'iYearwrk', 'dPOB', 'iSubfam2', 'iSex', 'iRownchld', 'iRiders', 'iDisabl2', 'dRearning', 'iYearsch', 'iTmpabsnt', 'iRemplpar', 'dIncome1', 'dOccup', 'dTravtime', 'dHours', 'iRlabor', 'dYrsserv', 'dRpincome', 'iMilitary', 'dWeek89', 'iWork89', 'dHour89', 'iPerscare', 'iWorklwk', 'iImmigr', 'iRrelchld', 'dIndustry', 'iMobillim'}

df.drop(labels=correlated_features, axis=1, inplace=True)
df.shape

(24583, 39)

Po obróbce zostało 38 kolumn w ramce. Jako, że wszystkie z nich są kategoryczne, zbadamy teraz ich rozkłady na grupowym histogramie.

Redukcja zmiennych: algorytm MCA

df1 = X.sample(frac=0.01, random_state=17)
df1.drop('caseid', inplace = True, axis = 1)

df_mca = mca.MCA(df1)

Sprawdźmy najpierw do jakiej liczby zmiennych możemy zredukować:

plt.plot(np.cumsum(df_mca.expl_var(greenacre=False, N=68)))
plt.xlabel('number of components')
plt.ylabel('cumulative explained variance');

Z wykresu wynika, że najbardziej odpowiednią wartością będzie 10.

mca = prince.MCA(n_components=10)
mca = mca.fit(df1)
mca = mca.transform(df1)
mca

	0	1	2	3	4	5	6	7	8	9
1551174	-0.695760	-0.366005	-0.118154	-0.232889	-0.092356	0.174044	0.000954	0.087259	-0.093411	-0.119080
1487409	1.077608	-0.013721	-0.039326	-0.032668	-0.014635	0.063703	0.053992	-0.033713	-0.036856	-0.030213
436156	-0.648278	0.923192	-0.427529	-0.059384	-0.139972	0.159919	-0.012562	0.067955	-0.045496	-0.023146
1088689	-0.695765	-0.460353	0.121820	-0.019130	-0.067644	0.107694	0.007207	0.031381	-0.134219	-0.017786
678665	-0.561437	0.819894	-0.151287	-0.156990	-0.084995	-0.053536	-0.403601	0.055249	0.300525	0.011119
...	...	...	...	...	...	...	...	...	...	...
77165	-0.651350	0.868671	-0.389745	-0.298434	-0.183215	0.337137	0.045245	0.046961	-0.005184	-0.022022
823271	-0.615616	0.668123	-0.125566	-0.235723	-0.191776	0.367406	0.012391	-0.011987	0.047139	-0.065916
1797933	-0.651221	0.496347	-0.106328	-0.112082	-0.045619	-0.654337	0.081591	0.146347	-0.834846	0.908348
1996241	-0.679699	-0.315659	0.067654	-0.178622	-0.117867	0.270970	-0.043009	0.055422	-0.072660	-0.099030
479824	-0.642057	-0.282944	-0.163976	0.558653	0.347452	-0.335598	0.242525	-0.247698	-0.033250	-0.432841

24583 rows × 10 columns

Redukcja zmiennych: współczynnik V Craméra

def cramers_v(x, y):
    confusion_matrix = pd.crosstab(x,y)
    chi2 = ss.chi2_contingency(confusion_matrix)[0]
    n = confusion_matrix.sum().sum()
    phi2 = chi2/n
    r,k = confusion_matrix.shape
    phi2corr = max(0, phi2-((k-1)*(r-1))/(n-1))
    rcorr = r-((r-1)**2)/(n-1)
    kcorr = k-((k-1)**2)/(n-1)
    return np.sqrt(phi2corr/min((kcorr-1),(rcorr-1)))

correlated_features = set()
for i in range(len(df1.columns)):
    for j in range(i):
        if cramers_v(df1[df1.columns[i]], df1[df1.columns[j]]) >= 0.70:
            colname = df1.columns[i]
            correlated_features.add(colname)
            
print(correlated_features)

{'iYearwrk', 'iRvetserv', 'dPOB', 'iSex', 'iWWII', 'iMobility', 'iRownchld', 'iVietnam', 'iDisabl2', 'dRearning', 'iYearsch', 'iTmpabsnt', 'iRspouse', 'dIncome5', 'dIncome2', 'iSept80', 'iRlabor', 'dYrsserv', 'dRpincome', 'iMilitary', 'dWeek89', 'iWork89', 'iPerscare', 'iLang1', 'iWorklwk', 'iRPOB', 'iSchool', 'iRrelchld', 'iMobillim'}

cols = df1.columns.values.tolist()
corrM = np.zeros((len(cols),len(cols)))
# there's probably a nice pandas way to do this
for col1, col2 in itertools.combinations(cols, 2):
    idx1, idx2 = cols.index(col1), cols.index(col2)
    corrM[idx1, idx2] = cramers_v(df1[col1], df1[col2])
    corrM[idx2, idx1] = corrM[idx1, idx2]

corr = pd.DataFrame(corrM, index=cols, columns=cols)
fig, ax = plt.subplots(figsize=(40, 40))
ax = sns.heatmap(corr, annot=True, ax=ax); ax.set_title("Cramer V Correlation between Variables");

powrót do EDA

df.hist(figsize=(15,20))
plt.show()

Możemy zauważyć, że większość zmiennych ma bardzo nierównomierny rozkład, w którym przeważa jedna lub dwie wartości. Jednynie zmienna dAge przypomina rozkład jednostajny.

Następnie przyjrzemy się wykresom wizualizującym relacje między zmiennymi, aby odkryć jakieś zależności, które mogą nam pomóc przy doborze klastrów. Przyjrzemy się bubble plotom, które pokazują liczności przecięć dwóch zmiennych.

def HeatMap(column1, column2):
    df1 = df.groupby([column1, column2]).size().reset_index(name='counts')
    df1 = pd.pivot(df1, index = 'race', columns = 'sex', values = 'counts')
    sns.heatmap(df1, annot=True, cmap='Reds', ax=ax1, fmt='g', norm=LogNorm())
    
def BubblePlot(column1, column2):
    df1 = df.groupby([column1, column2]).size().reset_index(name='counts')
    plt.figure(figsize = (12, 8))
    plt.scatter(df1[column1], df1[column2], s=df1['counts'], alpha = 0.7)
    plt.xlabel(column1)
    plt.gca().xaxis.set_major_locator(plt.MultipleLocator(1))
    plt.gca().yaxis.set_major_locator(plt.MultipleLocator(1))
    plt.ylabel(column2)

df1 = df.groupby(['iMarital', 'iClass']).size().reset_index(name='counts')

BubblePlot('iMarital', 'iClass')

Legenda:

	iMarital (marital status)	iClass (class of worker)
0	Now Married, Except Separated	N/a Less Than 16 Yrs. Old/unemp. Who Never Worked
1	Widowed	Emp. of a Private for Profit Company
2	Divorced	Emp. of a Private Not for Profit
3	Separated	Local Gov. Emp. City, County, Etc.
4	Never Married or Under 15 Yrs. Old	State Gov. Emp.
5		Federal Gov. Emp.
6		Self Emp. in Own Not Incorp.d Business
7		Self Emp. in Own Incorp.d Business
8		Working Without Pay in Family Business or Farm
9		Unemp., Last Worked in 1984 or Earlier

Widać, że przeważają grupy osób, które są w związku małżeńskim (1 na osi x) nigdy nie były w sformalizowanym związku lub poniiżej 15 roku życia (4 na osi x). Dodatkowo widać, że w grupie 4 na osi x, jest dużo osób, które nigdy nie były zatrudnione lub mają mniej niż 16 lat (0 na osi y), co jest dosyć intuicyjne. Dodatkowo wśród osób zamężnych jest przewaga osób pracujących prywatnie (1 na osi y)

BubblePlot('dAge', "iRelat1")

Legenda:

	iRelat1 (related or not)	dAge
0	Householder	0
1	Husband/wife	<13
2	Son/daughter	13-19
3	Stepson/stepdaughter	20-29
4	Brother/sister	30-39
5	Father/mother	40-49
6	Grandchild	50-64
7	Other Rel.	>65
8	Roomer/boarder/foster Child
9	Housemate/roommate
10	Unmarried Partner
11	Other Nonrel.
12	Instit. Person
13	Other Pers. in Grp. Qtrs.

Na wykresie ze zmiennymi dAge i iRelat1 widać niesamowicie ciekawe i oszałamiające rzeczy. Widać, że dorosłe osoby (oznaczone numerami większymi niż 2 na osi x) są głowami rodziny (0 na osi y) i są zakawalifikowani jako mąż lub żona (1 na osi y). Natomiast osoby młodsze (1 na osi x) są kwalifikowani przede wszytskim jako dzieci, co jest spójne.

BubblePlot('iEnglish', 'iSchool')

Legenda:

	iEnglish (ability to speak English)	iSchool (school enrollment)
0	N/a Less Than 5 Yrs. Old/speaks Only English	N/a Less Than 3 Yrs. Old
1	Very well	Not Attending School
2	Well	Yes, Pub. School, Pub. Coll.
3	Not Well	Yes, Private School, Private Coll.
4	Not at all

Największa grupą wśród rozważanych stanową osoby, które mówią tylko po angielsku. Wśród nich najwięcej jest osób, które nie chodzą do szkoły, czyli jak można wnioskować - dorosłych. Można też zauważyć, że więcej osób chodzi do publicznych, niż prywatnych placówek.

Skalowanie

Aby znormalizować naszą ramkę danych, postanowiliśmy przeskalować dane, aby zawierały się tylko w zbiorze [0,1]. Zmienne nadal pozostaną dyksretne, ale będzie zdecydowanie łatwiej analizować dobierane klastry w późniejszych krokach. Mamy na uwadze również to, że takie skalowanie zmieni dobór klastrów, bo odległości między punktami będą ewaluowane inaczej.

for col in df.columns:
    min_max_scaler1 = MinMaxScaler()
    scaled_col1 = min_max_scaler1.fit_transform(df[col].values.reshape(-1, 1))
    df[col] = scaled_col1

Modelowanie

from sklearn.cluster import KMeans, AgglomerativeClustering
from sklearn.metrics import silhouette_score, davies_bouldin_score
from scipy.cluster import hierarchy
from sklearn import cluster, datasets, mixture
from sklearn_extra.cluster import KMedoids
from yellowbrick.cluster import SilhouetteVisualizer

clus_num = range(1, 20)
scores = []

for k in clus_num:
    model = KMeans(n_clusters = k)
    model.fit(df)
    scores.append(model.score(df) * -1)
    
plt.plot(clus_num, scores, '-o')
plt.xlabel('number of clusters, k')
plt.ylabel('Within-cluster sum of squares')
plt.xticks(clus_num)
plt.show()

Na początek zbadaliśmy metodą łokcia, jaka liczba klastrów będzie optymalna. Wykres scorów dla liczby klastrów od 1 do 19 jest słabo stromy, więc trudno wybrać jednoznacznie punkt przegięcie. Ustalmy jednak, ze znajduje się on w okolicach wartośći od 5 d 9

Uznaliśmy, ze dobrą metodą na wizualizacje jak model dzieli dane na klastry, jest wykorzystanie centroidów i porównywanie ich w zależności od klastra dla każdej ze zmiennych.

def CentroidsPlot(centroids):
    cent_df = pd.DataFrame(columns= ['cluster', 'variable', 'centroid'])
    rows = centroids.shape[0]
    clusters = [x+1 for x in range(rows)]
    for colname in centroids.columns:
        data = {'cluster' : clusters,
               'variable' : [colname for x in range(rows)],
               'centroid' : centroids[colname]}
        cent_df = pd.concat([cent_df, pd.DataFrame(data)])
    plt.figure(figsize=(18, 8))
    scatter = plt.scatter(cent_df.variable, cent_df.centroid, c=cent_df.cluster, cmap='viridis', s = 200, alpha = 0.7)
    plt.legend(handles=scatter.legend_elements()[0], 
           labels=clusters,
           title="cluster")
    plt.xticks(rotation = 90, fontsize = 15)

KMeans

Na początek zbudujemy modele KMeans dla znalezionych wcześniej optymalnych liczb klastrów

model = KMeans(n_clusters = 5)
model.fit(df)
labs = model.labels_
centroids = model.cluster_centers_
centroids = pd.DataFrame(columns = df.columns, data = centroids)

CentroidsPlot(centroids)

model = KMeans(n_clusters = 7)
model.fit(df)
labs = model.labels_
centroids = model.cluster_centers_
centroids = pd.DataFrame(columns = df.columns, data = centroids)

CentroidsPlot(centroids)

model = KMeans(n_clusters = 9)
model.fit(df)
labs = model.labels_
centroids = model.cluster_centers_
centroids = pd.DataFrame(columns = df.columns, data = centroids)

CentroidsPlot(centroids)

Widać, że dobór większej liczby klastrów sprawia problemy, ponieważ kolejne klastry są do siebie za bardzo zbliżone. Dla 5 klastrów widać, że grupy mają różne dobrane centroidy, przez co wydaje się, że rzeczywiście dobrze rozdzielają rekory na różne klastry. Na pewno widać, że klaster oznaczony jako 3 dobrze oddziela młodzież i dzieci, natomiast klaster 4 ludzi starszych.